#include "geographicutil.h"
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/foreach.hpp>
#include <QtMath>
#include <QString>
#include <QVector>
#include <QDebug>
#include "shipAlgorithm/Mercator.h"
#include "geometricutil.h"
#include "commonutil.h"
#include "LonLatUtil.h"

using namespace ns_yht;

boostgm::model::polygon<GeoPoint> GeographicUtil::static_poly;
const double GeographicUtil::MetersPerNM = 1852.0;
const double GeographicUtil::MY_PI = 3.14159265359;

namespace
{
#define EARTH_RADIUS 6378137.0
#define EARTH_RADIUS2 6356752.3142451793

    typedef boostgm::model::point<double, 2, boostgm::cs::geographic<boostgm::degree>> GeoPoint;
    boostgm::srs::spheroid<double> spheroid(EARTH_RADIUS, EARTH_RADIUS2);

#define NORM(X)           \
    {                     \
        if (X < 0)        \
            X += 360;     \
        else if (X > 360) \
            X -= 360;     \
    }
}
GeographicUtil::GeographicUtil()
{
}

double GeographicUtil::getGeoArea(const DoublePairs &coordinates, bool isGeodesic)
{
    int n = coordinates.count();
    if (n < 3)
        return 0;
#if 0 
    bg::model::polygon<GeoPoint> geo_poly;
    for (int i = 0; i < coordinates.size(); i++)
    {
        GeoPoint p = {coordinates.at(i).second, coordinates.at(i).first};
        geo_poly.outer().push_back(p);
    }

    // Create geographic strategy with WGS84 spheroid
    // bg::srs::spheroid<double> spheroid(EARTH_RADIUS, EARTH_RADIUS2);
    // bg::strategy::area::geographic<> geo_strategy(spheroid);

    // Calculate the area of a geographic polygon
    // double area = qAbs(bg::area(geo_poly, geo_strategy));
    double area = qAbs(bg::area(geo_poly));

    return area;
#endif

    QString path = QString("POLYGON((%1 %2").arg(coordinates[0].second).arg(coordinates[0].first);
    for (int i = 1; i < coordinates.count(); i++)
    {
        path += QString(",%1 %2").arg(coordinates[i].second).arg(coordinates[i].first);
    }
    path += "))";

    // Create spherical polygon
    // bg::model::polygon<bg::model::point<double, 2, bg::cs::spherical_equatorial<bg::degree> > > sph_poly;

    // bg::read_wkt("POLYGON((123 30,122 28,124 29,122 29,124 28,123 30))", sph_poly);

    // Create spherical strategy with mean Earth radius in meters
    // bg::strategy::area::spherical<> sph_strategy(6371008.8);
    // bg::strategy::area::spherical<> sph_strategy(6377830);

    // Calculate the area of a spherical polygon
    //        double area = bg::area(sph_poly, sph_strategy);
    //        std::cout << "Area: " << area << std::endl;

    // Create geographic polygon
    boostgm::model::polygon<boostgm::model::point<double, 2, boostgm::cs::geographic<boostgm::degree>>> geo_poly;
    boostgm::read_wkt(path.toStdString(), geo_poly);

    // Create geographic strategy with WGS84 spheroid
    boostgm::strategy::area::geographic<> geo_strategy(spheroid);

    // Calculate the area of a geographic polygon
    double area = qAbs(boostgm::area(geo_poly, geo_strategy));

    return area;
}

double GeographicUtil::getGeoDistance(const QVector<QPair<double, double>> &coordinates, bool isGeodesic)
{
    int n = coordinates.count();
    if (n < 2)
        return 0;

    // Create geographic strategy with WGS84 spheroid
    // bg::srs::spheroid<double> spheroid(EARTH_RADIUS, EARTH_RADIUS2);
    // bg::strategy::distance::geographic<> geo_strategy(spheroid);

    GeoPoint pt1(coordinates[0].second, coordinates[0].first);
    double distance = 0;
    for (int i = 1; i < coordinates.count(); i++)
    {
        GeoPoint pt2(coordinates[i].second, coordinates[i].first);
        // distance += bg::distance(pt1, pt2, geo_strategy);
        distance += boostgm::distance(pt1, pt2);
        pt1 = pt2;
    }

    return distance;
#if 0
    double distance = 0;
    for (int i = 1; i < coordinates.count(); i++)
    {
        double s12, azi;
        if (isGeodesic)
            geod.Inverse(coordinates.at(i - 1).first, coordinates.at(i - 1).second, coordinates.at(i).first, coordinates.at(i).second, s12);
        else
            rhumb.Inverse(coordinates.at(i - 1).first, coordinates.at(i - 1).second, coordinates.at(i).first, coordinates.at(i).second, s12, azi);
        distance += s12;
    }

    return distance;
#endif

#if 0
    double area = 0.0;
    double perimeter = 0.0;
    if (isGeodesic)
    {
        PolygonArea polygon(geod);
        for (int i = 0; i < n; i++)
            polygon.AddPoint(coordinates.at(i).first, coordinates.at(i).second);

        polygon.Compute(false, true, perimeter, area);
    }
    else
    {
        PolygonAreaT<Rhumb> polygon(rhumb);
        for (int i = 0; i < n; i++)
            polygon.AddPoint(coordinates.at(i).first, coordinates.at(i).second);
        polygon.Compute(false, true, perimeter, area);
    }
    return perimeter;   //这里计算的是周长，如果只有两个点，相当于是两倍
#endif
}

double GeographicUtil::getGeoDistance(const DoublePair &point, const DoublePair &segPoint0, const DoublePair &segPoint1)
{
    using GeoSegment = boostgm::model::segment<GeoPoint>;

    GeoPoint pt = {point.second, point.first};
    GeoSegment seg = {{segPoint0.second, segPoint0.first}, {segPoint1.second, segPoint1.first}};

    return boostgm::distance(pt, seg);
}

double GeographicUtil::getGeoDistance(const QPair<double, double> &coor1, const QPair<double, double> &coor2, bool isGeodesic)
{

    GeoPoint pt1(coor1.second, coor1.first);
    GeoPoint pt2(coor2.second, coor2.first);
    return boostgm::distance(pt1, pt2);
#if 0
    if (isGeodesic)
        return calculateGeodesicDistance(coor1, coor2) * MetersPerNM;
    else
        return calculateRhumbDistance(coor1, coor2);
#endif
}
/**
 * @brief 计算一个点到一个多边形的距离
 *
 * @param point 点坐标
 * @param shape  多边形
 * @return double  单位：米
 */
double GeographicUtil::distanceToPolygon(const DoublePair &point, const DoublePairs &shape)
{
    GeoPoint pt = {point.second, point.first};

    boostgm::model::polygon<GeoPoint> geo_poly;
    for (int i = 0; i < shape.count(); i++)
    {
        GeoPoint p = {shape[i].second, shape[i].first};
        geo_poly.outer().push_back(p);
    }
    return boostgm::distance(pt, geo_poly);
}

double GeographicUtil::distanceToPolyline(const DoublePair &point, const DoublePairs &shape)
{
    GeoPoint pt = {point.second, point.first};
    boostgm::model::linestring<GeoPoint> geo_line;
    for (int i = 0; i < shape.count(); i++)
    {
        GeoPoint p = {shape[i].second, shape[i].first};
        geo_line.push_back(p);
    }
    return boostgm::distance(pt, geo_line);
}

double GeographicUtil::distanceToPie(const DoublePair &point, const DoublePairs &shape)
{
    GeoPoint pt = {point.second, point.first};
    boostgm::model::ring<GeoPoint> geo_ring;
    for (int i = 0; i < shape.count(); i++)
    {
        GeoPoint p = {shape[i].second, shape[i].first};
        geo_ring.push_back(p);
    }
    return boostgm::distance(pt, geo_ring);
}

double GeographicUtil::getGeoAngle(const QPair<double, double> &coor1, const QPair<double, double> &coor2)
{
    GeoPoint p1(coor1.second, coor1.first);
    GeoPoint p2(coor2.second, coor2.first);
    return qRadiansToDegrees(boostgm::azimuth(p1, p2));

#if 0
    double s01, azi;
    rhumb.Inverse(coor1.first, coor1.second, coor2.first, coor2.second, s01, azi);
    return azi;
#endif
}

/**
    这两个方法应该是等价的。
    如果你的需求是计算两个经纬度点之间的方位角度，且精度要求不高，那么使用三角函数和反正切函数的方法足够准确，并且计算较为简单快速。
    然而，如果你的需求是计算两个经纬度点在地球表面上的实际方位角度，且对精度有较高要求，那么转换为墨卡托坐标再计算方位角的方法可能更准确,
    也就是上面的方法.
    转换为墨卡托坐标可以将地球表面上的点映射到一个平面上，从而消除了地球曲面引起的误差，并且可以使用更精确的数学模型来计算方位角度。
    总之，如果你的应用场景对精度要求不高，或者只是需要一个近似的方位角度，那么使用三角函数和反正切函数的方法是足够准确的。
    如果你的应用场景对精度要求较高，或者需要在地球表面上进行准确的方位角度计算，那么转换为墨卡托坐标再计算方位角的方法可能更适合。
**/
#if 0
double GeographicUtil::getGeoAngle(const QPair<double, double> &coor1, const QPair<double, double> &coor2)
{
    GeoPoint p1(coor1.second, coor1.first);
    GeoPoint p2(coor2.second, coor2.first);

    return qRadiansToDegrees(bg::azimuth(p1, p2));
}


double calculateBearing(double lat1, double lon1, double lat2, double lon2)
{
    double dLon = toRadians(lon2 - lon1);
    double y = sin(dLon) * cos(toRadians(lat2));
    double x = cos(toRadians(lat1)) * sin(toRadians(lat2)) - sin(toRadians(lat1)) * cos(toRadians(lat2)) * cos(dLon);
    double bearing = atan2(y, x);
    bearing = fmod((bearing + 2 * M_PI), (2 * M_PI)); // Convert bearing to a positive value between 0 and 2*pi
    bearing = bearing * 180.0 / M_PI;                 // Convert bearing to degrees
    return bearing;
}
#endif

double GeographicUtil::getGeoAngle_Radians(const DoublePair &coor1, const DoublePair &coor2)
{
    GeoPoint p1(coor1.second, coor1.first);
    GeoPoint p2(coor2.second, coor2.first);
    return boostgm::azimuth(p1, p2);
}

double GeographicUtil::getLat(double lat1, double distance, bool up)
{
    if (up)
        return qRadiansToDegrees(qDegreesToRadians(lat1) + distance / EARTH_RADIUS);
    else
        return qRadiansToDegrees(qDegreesToRadians(lat1) - distance / EARTH_RADIUS);
}

double GeographicUtil::getLon(double lat1, double lon1, double distance, bool up)
{

    using segment_type = boostgm::model::segment<GeoPoint>;

    boostgm::strategy::line_interpolate::geographic<boostgm::strategy::vincenty> geo_strategy(spheroid);

    GeoPoint pt1(lon1, lat1);
    GeoPoint pt2(10 + lon1, lat1);
    segment_type const s{pt1, pt2};

    GeoPoint p;
    line_interpolate(s, distance, p, geo_strategy);

    double lon2 = p.get<0>();
    if (!up)
        lon2 = lon1 * 2 - lon2;
    return lon2;

#if 0
    RhumbLine line = rhumb.Line(lat1, lon1, up ? 90 : -90);
    double lat2, lon2;
    line.Position(distance, lat2, lon2);
    return lon2;
#endif
}

void GeographicUtil::getCoordinate(double lat1, double lon1, double angle, double distance, double &lat2, double &lon2)
{
    if (angle < 0)
        angle += 360;
    else if (angle > 360)
        angle -= 360;

    double radian = qDegreesToRadians(angle);

    double disx = qAbs(distance * qSin(radian));
    bool upx = (angle >= 0 && angle < 180);
    lon2 = getLon(lat1, lon1, disx, upx);

    double disy = qAbs(distance * qCos(radian));
    bool upy = !(angle > 90 && angle < 270);
    lat2 = getLat(lat1, disy, upy);

#if 0
    RhumbLine line = rhumb.Line(lat1, lon1, angle);
    line.Position(distance, lat2, lon2);
#endif
}

void GeographicUtil::getCoordinate(const QPair<double, double> &coor1, const QPair<double, double> &coor2, double distance, QPair<double, double> &newCoor)
{

    double dis1 = getGeoDistance(coor1, coor2);
    double scale = distance / dis1;
    newCoor.second = (coor2.second - coor1.second) * scale + coor1.second;
    newCoor.first = (coor2.first - coor1.first) * scale + coor1.first;

#if 0
    double s12, azi;
    rhumb.Inverse(coor1.first, coor1.second, coor2.first, coor2.second, s12, azi);
    RhumbLine line = rhumb.Line(coor1.first, coor1.second, azi);
    line.Position(distance, newCoor.first, newCoor.second);
#endif
}

void GeographicUtil::simplifyLine(const QVector<QPair<double, double>> &spoints, QVector<QPair<double, double>> &dpoints, double dis)
{
    typedef boost::geometry::model::d2::point_xy<double> xy;

    boost::geometry::model::linestring<xy> line;
    for (int i = 0; i < spoints.count(); i++)
        line.push_back(xy(spoints[i].second, spoints[i].first));

    // Simplify it, using distance of 0.5 units
    boost::geometry::model::linestring<xy> simplified;
    boost::geometry::simplify(line, simplified, dis);

    for (boost::geometry::model::linestring<xy>::iterator it = simplified.begin(); it != simplified.end(); it++)
    {
        dpoints.append(QPair<double, double>(it->y(), it->x()));
    }
}

void GeographicUtil::setStaticPolygon(const DoublePairs &coordinates)
{
    static_poly.clear();
    int n = coordinates.count();
    if (n < 3)
        return;
    for (int i = 0; i < n; i++)
    {
        GeoPoint p(coordinates[i].second, coordinates[i].first);
        static_poly.outer().push_back(p);
    }
    if (qAbs(coordinates.first().first - coordinates.last().first) > 1e-6 || qAbs(coordinates.first().second - coordinates.last().second) > 1e-6)
        static_poly.outer().push_back(static_poly.outer().front());
}

bool GeographicUtil::pointInGeometry(double lat, double lon)
{
    return boostgm::within(GeoPoint(lon, lat), static_poly);
}

bool GeographicUtil::pointInGeometry(double lat, double lon, const QVector<QPair<double, double>> &coordinates)
{
    int n = coordinates.count();
    if (n < 3)
        return 0;

#if 0
    QString path = QString("POLYGON((%1 %2").arg(coordinates[0].second).arg(coordinates[0].first);
    for (int i = 1; i < coordinates.count(); i++)
    {
        path += QString(",%1 %2").arg(coordinates[i].second).arg(coordinates[i].first);
    }
    path += QString(",%1 %2").arg(coordinates[0].second).arg(coordinates[0].first);
    path += "))";
    return pointInGeometry(lat, lon, path);
#endif

    boostgm::model::polygon<GeoPoint> geo_poly;
    for (int i = 0; i < n; i++)
    {
        GeoPoint p(coordinates[i].second, coordinates[i].first);
        geo_poly.outer().push_back(p);
    }
    geo_poly.outer().push_back(geo_poly.outer().front());
    return boostgm::within(GeoPoint(lon, lat), geo_poly);
}

bool GeographicUtil::pointInGeometry(double lat, double lon, const QString &coordinates)
{
    boostgm::model::polygon<GeoPoint> geo_poly;
    boostgm::read_wkt(coordinates.toStdString(), geo_poly);

    return boostgm::within(GeoPoint(lon, lat), geo_poly);
}

bool GeographicUtil::intersects_lineToPoly(const DoublePairs &coordinates1)
{
    boostgm::model::linestring<GeoPoint> geo_poly1;
    for (int i = 0; i < coordinates1.count(); i++)
    {
        GeoPoint p(coordinates1.at(i).second, coordinates1.at(i).first);
        geo_poly1.push_back(p);
    }
    return boostgm::intersects(geo_poly1, static_poly);
}

bool GeographicUtil::intersects_lineToPoly(const DoublePairs &coordinates1, const DoublePairs &coordinates2)
{
    boostgm::model::linestring<GeoPoint> geo_poly1;
    for (int i = 0; i < coordinates1.count(); i++)
    {
        GeoPoint p(coordinates1.at(i).second, coordinates1.at(i).first);
        geo_poly1.push_back(p);
    }

    boostgm::model::polygon<GeoPoint> geo_poly2;
    for (int i = 0; i < coordinates2.count(); i++)
    {
        GeoPoint p(coordinates2.at(i).second, coordinates2.at(i).first);
        geo_poly2.outer().push_back(p);
    }
    return boostgm::intersects(geo_poly1, geo_poly2);
}

bool GeographicUtil::intersects_polyToPoly(const DoublePairs &coordinates1)
{
    boostgm::model::polygon<GeoPoint> geo_poly1;
    for (int i = 0; i < coordinates1.count(); i++)
    {
        GeoPoint p(coordinates1.at(i).second, coordinates1.at(i).first);
        geo_poly1.outer().push_back(p);
    }
    return boostgm::intersects(geo_poly1, static_poly);
}

bool GeographicUtil::intersects_polyToPoly(const DoublePairs &coordinates1, const DoublePairs &coordinates2)
{
    boostgm::model::polygon<GeoPoint> geo_poly1;
    for (int i = 0; i < coordinates1.count(); i++)
    {
        GeoPoint p(coordinates1.at(i).second, coordinates1.at(i).first);
        geo_poly1.outer().push_back(p);
    }

    boostgm::model::polygon<GeoPoint> geo_poly2;
    for (int i = 0; i < coordinates2.count(); i++)
    {
        GeoPoint p(coordinates2.at(i).second, coordinates2.at(i).first);
        geo_poly2.outer().push_back(p);
    }
    return boostgm::intersects(geo_poly1, geo_poly2);
}

bool GeographicUtil::intersection_polyToPoly(const QVector<DoublePairs> &geometry1, const QVector<DoublePairs> &geometry2, QVector<DoublePairs> &result)
{
    result.clear();
    boostgm::model::polygon<GeoPoint> geo_poly1;
    convertPolygonToPolygon<boostgm::model::polygon<GeoPoint>>(geometry1, geo_poly1);

    boostgm::model::polygon<GeoPoint> geo_poly2;
    convertPolygonToPolygon<boostgm::model::polygon<GeoPoint>>(geometry2, geo_poly2);
    std::deque<boostgm::model::polygon<GeoPoint>> output;
    bool res = boostgm::intersection(geo_poly1, geo_poly2, output);
    if (output.empty())
        return false;
    BOOST_FOREACH (boostgm::model::polygon<GeoPoint> const &p, output)
    {
        DoublePairs dpoints;
        for (int i = 0; i < p.outer().size(); i++)
        {
            dpoints.append(QPair<double, double>(p.outer()[i].get<1>(), p.outer()[i].get<0>()));
        }
        result.append(dpoints);
    }
    return true;
}

template <class T>
void GeographicUtil::convertPolygonToPolygon(const QVector<DoublePairs> &geometry, T &result)
{
    result.clear();
    QString path = QString("POLYGON(");
    for (int t = 0; t < geometry.count(); t++)
    {
        auto coors = geometry[t];
        int n = coors.count();
        if (n == 0)
            continue;
        QString s = QString("(%1 %2").arg(coors[0].second, 0, 'f', 10).arg(coors[0].first, 0, 'f', 10);

        for (int i = 1; i < n; i++)
        {
            s += QString(",%1 %2").arg(coors[i].second, 0, 'f', 10).arg(coors[i].first, 0, 'f', 10);
        }
        s += ")";
        path += s;
    }
    path.append(")");
    boostgm::read_wkt(path.toStdString(), result);
}

void GeographicUtil::getPointCircle(double lat, double lon, double distance, DoublePairs &circle, int circlePointCount)
{
    // Declare the geographic_point_circle strategy (with 36 points)
    // Default template arguments (taking Andoyer strategy)
    boostgm::strategy::buffer::geographic_point_circle<> point_strategy(circlePointCount);

    // Declare the distance strategy (one kilometer, around the point, on Earth)
    boostgm::strategy::buffer::distance_symmetric<double> distance_strategy(distance);

    // Declare other necessary strategies, unused for point
    boostgm::strategy::buffer::join_round join_strategy;
    boostgm::strategy::buffer::end_round end_strategy;
    boostgm::strategy::buffer::side_straight side_strategy;

    // Declare/fill a point on Earth, near Amsterdam
    GeoPoint p(lon, lat);

    // Create the buffer of a point on the Earth
    boostgm::model::multi_polygon<boostgm::model::polygon<GeoPoint>> result;
    boostgm::buffer(p, result,
                    distance_strategy, side_strategy,
                    join_strategy, end_strategy, point_strategy);

    for (size_t i = 0, n = result.at(0).outer().size(); i < n; i++)
    {
        circle.append(DoublePair(result.at(0).outer().at(i).get<1>(), result.at(0).outer().at(i).get<0>()));
    }
}

void GeographicUtil::getParallelLine(const DoublePairs &line, double distance, bool isLeft, DoublePairs &parallel)
{
    QVector<QPointF> lineXY;
    for (int i = 0; i < line.count(); i++)
    {
        lineXY.append(QPointF(Mercator::ToX(line[i].second), Mercator::ToY(line[i].first)));
    }

    QVector<QPointF> outLine;
    GeometricUtil::parallelLine(lineXY, distance, isLeft, outLine);
    parallel.clear();
    for (int i = 0; i < outLine.count(); i++)
    {
        parallel.append(DoublePair(Mercator::ToLat(outLine.at(i).y()), Mercator::ToLon(outLine.at(i).x())));
    }
}

void GeographicUtil::getParallelLines(const DoublePairs &line, double loffset, double roffset, DoublePairs &lParallel, DoublePairs &rParallel)
{
    getParallelLine(line, loffset, true, lParallel);
    getParallelLine(line, roffset, false, rParallel);
}

void GeographicUtil::getArc(double lat, double lon, double distance, double startAngle, double sweepAngle, DoublePairs &arc)
{

    int arcPointCount = 36;
    // Calculate the points on a circle centered at (lat, lon) with radius distance
    DoublePairs arc1;
    getPointCircle(lat, lon, distance, arc1, arcPointCount);

    // 获取起始点
    DoublePair center(lat, lon);
    double minDelta = 99999;
    int minI = 0;
    for (int i = 0; i < arc1.size(); i++)
    {
        double a = getGeoAngle(center, arc1.at(i));
        double delta = qAbs(a - startAngle);
        if (delta < minDelta)
        {
            minDelta = delta;
            minI = i;
        }
    }

    // Generate the arc points based on the startAngle and sweepAngle
    arc.clear();
    double sw = 0;
    double step = 360.0 / arcPointCount;
    bool finished = false;
    for (int i = minI; i < arc1.size(); i++)
    {
        arc.append(arc1.at(i));
        sw += step;
        if (sw > sweepAngle)
        {
            finished = true;
            break;
        }
    }
    if (!finished)
    {
        for (int i = 0; i < minI; i++)
        {
            arc.append(arc1.at(i));
            sw += step;
            if (sw > sweepAngle)
            {
                break;
            }
        }
    }

#if 0
    // arc.clear();
    double lat2, lon2;
    if (sweepAngle > 0)
    {
        int step = 1;

        for (int i = 0; i < sweepAngle; i += step)
        {
            double angle = startAngle + i;
            rhumb.Direct(lat, lon, angle, distance, lat2, lon2);
            arc.append(DoublePair(lat2, lon2));
        }
    }
    else
    {
        int step = -1;
        for (int i = 0; i > sweepAngle; i += step)
        {
            double angle = startAngle + i;
            rhumb.Direct(lat, lon, angle, distance, lat2, lon2);
            arc.append(DoublePair(lat2, lon2));
        }
    }
    rhumb.Direct(lat, lon, startAngle + sweepAngle, distance, lat2, lon2);
    arc.append(DoublePair(lat2, lon2));
#endif
}

double GeographicUtil::calculateRhumbDistance(double lat1, double lon1, double lat2, double lon2)
{

    // 将经纬度转换为弧度
    lat1 = qDegreesToRadians(lat1);
    lon1 = qDegreesToRadians(lon1);
    lat2 = qDegreesToRadians(lat2);
    lon2 = qDegreesToRadians(lon2);

    // // Haversine 公式计算恒向线距离
    double dLat = abs(lat2 - lat1);
    double dLon = abs(lon2 - lon1);
    if (dLon > M_PI)
    {
        dLon = 2 * M_PI - dLon;
    }

    // 来自书本和网络资料
    // 同一纬度线上
    if (CommonUtil::floatEqual(lat1, lat2, 1e-6))
        return qRadiansToDegrees(abs(dLon * cos(lat1))) * 60;

    double t2 = tan(M_PI_4 + lat2 / 2);
    double t1 = tan(M_PI_4 + lat1 / 2);
    double c = atan2(dLon, log(t2) - log(t1));

    double s = abs(dLat / cos(c));
    return qRadiansToDegrees(s) * 60;

#if 0
    double s12, azi;
    rhumb.Inverse(lat1, lon1, lat2, lon2, s12, azi);
    return s12;
#endif
}

double GeographicUtil::calculateRhumbDistance(const DoublePair &point1, const DoublePair &point2)
{
    return calculateRhumbDistance(point1.first, point1.second, point2.first, point2.second);
}

double GeographicUtil::calculateHeading(double lat1, double lon1, double lat2, double lon2)
{

    lat1 = qDegreesToRadians(lat1);
    lon1 = qDegreesToRadians(lon1);
    lat2 = qDegreesToRadians(lat2);
    lon2 = qDegreesToRadians(lon2);
    double deltaLong = lon2 - lon1;
    double dLat = log((tan(lat2 / 2) + cos(lat1)) / (tan(lat1 / 2) + cos(lat2)));
    double dLon = deltaLong / 2 * sin(lat1 - lat2);
    double heading = qRadiansToDegrees(atan2(dLon, dLat));

    if (heading < 0)
        heading += 360.0;
    return heading;
#if 0
    double s12, azi;
    rhumb.Inverse(lat1, lon1, lat2, lon2, s12, azi);
    return azi;
#endif
}

double GeographicUtil::calculateHeading(const DoublePair &point1, const DoublePair &point2)
{
    return calculateHeading(point1.first, point1.second, point2.first, point2.second);
}

double GeographicUtil::calculateGeodesicDistance(double lat1, double lon1, double lat2, double lon2)
{
    // 将经纬度转换为弧度
    lat1 = qDegreesToRadians(lat1);
    lon1 = qDegreesToRadians(lon1);
    lat2 = qDegreesToRadians(lat2);
    lon2 = qDegreesToRadians(lon2);
    double dLon = abs(lon2 - lon1);
    if (dLon > M_PI)
    {
        dLon = 2 * M_PI - dLon;
    }
    // Haversine哈弗辛公式
    double cs = sin(lat1) * sin(lat2) + cos(lat1) * cos(lat2) * cos(dLon);
    double s = (acos(cs));
    double a = asin(cos(lat2) * sin(dLon) / sin(s));
    return qRadiansToDegrees(s) * 60; // 返回海里

    // double s12;
    // return geod.Inverse(lat1, lon1, lat2, lon2, s12) * 60;
}

double GeographicUtil::calculateGeodesicDistance(const DoublePair &point1, const DoublePair &point2)
{
    return calculateGeodesicDistance(point1.first, point1.second, point2.first, point2.second);
}

// void GeographicUtil::getEllipseCoordinates(double slat, double slon, double elat, double elon, Coordinates &coordinates, double alpha) const
// {
//     if (!LonLatUtil::isValid(slat, slon) || !LonLatUtil::isValid(elat, elon))
//         return;
//     QPointF spoint = pixelForCoordinate(Coordinate(slat, slon));
//     QPointF epoint = pixelForCoordinate(Coordinate(elat, elon));
//     QPointF basePoint = pixelForCoordinate(Coordinate((slat + elat) / 2, (slon + elon) / 2));

//     double sy = spoint.y();
//     double sx = spoint.x();
//     double ex = epoint.x();
//     double ey = epoint.y();

//     if (sx > ex)
//     {
//         double tmp = sx;
//         sx = ex;
//         ex = tmp;
//     }
//     if (sy > ey)
//     {
//         double tmp = sy;
//         sy = ey;
//         ey = tmp;
//     }

//     double width = sx - ex;
//     double height = sy - ey;

//     double r1 = width / 2, r2 = height / 2;
//     double cx = ex + r1, cy = ey + r2;

//     double a = sqrt(sqrt(width * width + height * height));
//     if (a < 1)
//         a = 1.0;
//     double steps = floor(a * M_PI);

//     for (double d = 0; d < 2 * M_PI - M_PI / steps; d += 2 * M_PI / steps)
//     {
//         QPointF point = QPointF(cx + r1 * sin(d), cy - r2 * cos(d));
//         point = GeometricUtil::rotateOnePoint(basePoint, alpha, point);
//         Coordinate coor = coordinateForPixel(point);
//         coordinates.push_back(coor);
//     }
//     coordinates.push_back(coordinates.first());
// }

void GeographicUtil::getArcCoordinates(const DoublePair &center, const DoublePair &startPoint, double sweepAngle, DoublePairs &coordinates)
{
    // if (!LonLatUtil::isValid(startPoint.first, startPoint.second) || !LonLatUtil::isValid(center.first, center.second))
    //     return;
    // double s12, azi;
    // GeographicUtil::rhumb.Inverse(center.first, center.second, startPoint.first, startPoint.second, s12, azi);
    // GeographicUtil::getArc(center.first, center.second, s12, azi, sweepAngle, coordinates);
}

void GeographicUtil::getCircleCoordinates(const DoublePair &center, const DoublePair &radiusPoint, DoublePairs &coordinates)
{
    if (!LonLatUtil::isValid(radiusPoint.first, radiusPoint.second) || !LonLatUtil::isValid(center.first, center.second))
        return;

    double distance = getGeoDistance(center, radiusPoint, false);
    getPointCircle(center.first, center.second, distance, coordinates, 36 * 5);

    // QPointF spoint = pixelForCoordinate(radiusPoint);
    // QPointF cpoint = pixelForCoordinate(center);

    // double radius = qSqrt(qPow((spoint.x() - cpoint.x()), 2.0) + qPow((spoint.y() - cpoint.y()), 2.0));

    // double cx = cpoint.x(), cy = cpoint.y();

    // double a = sqrt(sqrt(2 * radius * radius));
    // if (a < 1)
    //     a = 1.0;
    // double steps = floor(a * M_PI) * 10;

    // double startAngle = 0;
    // double endAngle = 2 * M_PI;
    // for (double a = startAngle; a < endAngle; a += 2 * M_PI / steps)
    // {
    //     QPointF point = QPointF(cx + radius * sin(a), cy - radius * cos(a));

    //     Coordinate coor = coordinateForPixel(point);
    //     coordinates.push_back(coor);
    // }
    // coordinates.push_back(coordinates.first());
}

void GeographicUtil::getCircleCoordinates(const DoublePair &center, double radius, DoublePairs &coordinates)
{
    getPointCircle(center.first, center.second, radius, coordinates, 36 * 5);
}

void GeographicUtil::getGreateCircleCoors(const DoublePairs &lineCoors, DoublePairs &greatCircleCoors, DoublePairs &centerCoors)
{
    // if (lineCoors.count() < 2)
    //     return;
    // Coordinate start = lineCoors[0];
    // Coordinate end = lineCoors[1];
    // Coordinates coorsseg;
    // GreatCircle::GetArc(start, end, greatCircleCoors);
    // if (greatCircleCoors.count() < 2)
    //     return;
    // centerCoors.append(greatCircleCoors.at(greatCircleCoors.count() / 2));
    // for (int i = 1; i < lineCoors.count() - 1; i++)
    // {
    //     start = lineCoors[i];
    //     end = lineCoors[i + 1];
    //     GreatCircle::GetArc(start, end, coorsseg);
    //     if (coorsseg.count() < 2)
    //         continue;
    //     centerCoors.append(coorsseg.at(coorsseg.count() / 2));
    //     coorsseg.pop_front();
    //     greatCircleCoors.append(coorsseg);
    // }

    if (lineCoors.count() < 2)
        return;
    // DoublePair start = lineCoors[0];
    // DoublePair end = lineCoors[1];

    // const double disStep = 1852 * 2;
    // greatCircleCoors.clear();
    // GeographicLib::GeodesicLine line = GeographicUtil::geod.InverseLine(start.first, start.second, end.first, end.second);
    // double s12 = line.Distance();
    // int n = qFloor(s12 / disStep);
    // for (int i = 1; i < n; i++)
    // {
    //     double s = i * disStep;
    //     double lat, lon, azi2;
    //     line.Position(s, lat, lon, azi2);
    //     greatCircleCoors.append({lat, lon});
    // }
    greatCircleCoors.append(lineCoors.at(0));
    for (int i = 1; i < lineCoors.count(); i++)
    {
        const DoublePair &A0 = lineCoors.at(i - 1);
        const DoublePair &A1 = lineCoors.at(i);

        // double s01, azi0, azi1;

        // geod.Inverse(A0.first, A0.second, A1.first, A1.second, s01, azi0, azi1);
        // NORM(azi0)
        // NORM(azi1)

        DoublePairs coors;
        getGreateCircleCoors(A0.first, A0.second, A1.first, A1.second, coors);
        greatCircleCoors.append(coors);
        greatCircleCoors.append(A1);
        if (coors.count() == 0)
            centerCoors.append(DoublePair((A1.first + A0.first) / 2, (A1.second + A0.second) / 2));
        else
            centerCoors.append(coors.at(coors.count() / 2));
    }
}

void GeographicUtil::getGreateCircleCoors(double lat1, double lon1, double lat2, double lon2, DoublePairs &greatCircleCoors)
{
#if 0
    const double disStep = 1852 * 2;
    greatCircleCoors.clear();
    GeodesicLine line = geod.InverseLine(lat1, lon1, lat2, lon2);
    double s12 = line.Distance();
    int n = qFloor(s12 / disStep);
    for (int i = 1; i < n; i++)
    {
        double s = i * disStep;
        double lat, lon, azi2;
        line.Position(s, lat, lon, azi2);
        greatCircleCoors.append({lat, lon});
    }
#endif
}

void GeographicUtil::intersects_line(const DoublePairs &coordinates1, const DoublePairs &coordinates2, DoublePair &intersect)
{
    GeoLineString geo_poly1;
    for (int i = 0; i < coordinates1.count(); i++)
    {
        GeoPoint p(coordinates1.at(i).second, coordinates1.at(i).first);
        geo_poly1.push_back(p);
    }

    GeoLineString geo_poly2;
    for (int i = 0; i < coordinates2.count(); i++)
    {
        GeoPoint p(coordinates2.at(i).second, coordinates2.at(i).first);
        geo_poly2.push_back(p);
    }
    std::vector<GeoPoint> geo;
    boostgm::intersection(geo_poly1, geo_poly2, geo);
    if (geo.size() < 1)
        return;
    for (int i = 0; i < geo.size(); i++)
    {
        qInfo() << "i" << i << geo[i].get<0>() << geo[i].get<1>();
    }
    intersect = (DoublePair(geo[0].get<1>(), geo[0].get<0>()));
}
